function [ret1, ret2] = LFD_C(image, radius, neighbors, Q, rng)
% Q is the number of disks, default is 20. It may depends on image size.
% RNG is the magnitude used, default is [2:neighbors/2+1].
%% preprocess
if nargin < 4
    Q = 20;
end
if nargin < 5
    rng = [2: neighbors / 2+1];
end

ret1 = [];
ret2 = [];

% grayscale
if size(image,3) ~= 1
    image = rgb2gray(image);
end
% normalize image
image = double(image);
image = image - mean(image(:)); %zero mean
image = image / var(image(:)) ^ 0.5; %unit variance

%% get pixel values for LBP
% code extracted from lbp
m = 1:neighbors;
spoints = [-radius*sin(2*pi*m/neighbors);...
    radius*cos(2*pi*m/neighbors)]'; % circular

miny=min(spoints(:,1));
maxy=max(spoints(:,1));
minx=min(spoints(:,2));
maxx=max(spoints(:,2));

% Block size, each LBP code is computed within a block of size bsizey*bsizex
bsizey=ceil(max(maxy,0))-floor(min(miny,0))+1;
bsizex=ceil(max(maxx,0))-floor(min(minx,0))+1;

% Coordinates of origin (0,0) in the block
origy=1-floor(min(miny,0));
origx=1-floor(min(minx,0));

[ysize xsize] = size(image);

% Calculate dx and dy;
dx = xsize - bsizex;
dy = ysize - bsizey;
d_image=image;

for i = 1:neighbors
    y = spoints(i,1)+origy;
    x = spoints(i,2)+origx;
    % Calculate floors, ceils and rounds for the x and y.
    fy = floor(y); cy = ceil(y); ry = round(y);
    fx = floor(x); cx = ceil(x); rx = round(x);
    % Check if interpolation is needed.
    if (abs(x - rx) < 1e-6) && (abs(y - ry) < 1e-6)
        % Interpolation is not needed, use original datatypes
        data(:, :, i) = image(ry:ry+dy,rx:rx+dx);
    else
        % Interpolation needed, use double type images 
        ty = y - fy;
        tx = x - fx;

        % Calculate the interpolation weights.
        w1 = (1 - tx) * (1 - ty);
        w2 =      tx  * (1 - ty);
        w3 = (1 - tx) *      ty ;
        w4 =      tx  *      ty ;
        % Compute interpolated pixel values
        data(:, :, i) = w1*d_image(fy:fy+dy,fx:fx+dx) + w2*d_image(fy:fy+dy,cx:cx+dx) + ...
        w3*d_image(cy:cy+dy,fx:fx+dx) + w4*d_image(cy:cy+dy,cx:cx+dx);
    end
end

%% FFT each set of points
[K,L,~] = size(data);
m = fft(data,[],3);
% only 1..N/2+1 are different
m = m(:, :, 1:floor(neighbors / 2) + 1);
% only use magitude
m = abs(m);

%% FFT2 each magnitude image
% calculate distance to center
centerK = (K+1)/2;
centerL = (L+1)/2;
dk = (repmat([1:K]', [1 L]) - centerK).^2;
dl = (repmat(1:L, [K 1]) - centerL).^2;
% d[i,j] = k only when distance in [k, k+1)
d = floor((dk + dl).^0.5)+1;
maxd = max(d(:));
t = histc(d(:),[1:maxd])'; % occurance    
lrng = length(rng);
for k = 1:lrng % skip first magnitude?why?
    i = rng(k);
    % FFT2
    p = m(:, :, i);
    p = fft2(p);
    p = fftshift(p);
    p = abs(p);
    
    % COUNT ALL DISKS, IN RADIUS (d,d+1) d in 0..maxd-1
    v = zeros(1, maxd);
    for j = 1:maxd
        u = p(d == j);
        v(j) = sum(u(:));
    end
        
    % special disks as features in paper, size = Q
    q = zeros(1, Q);
    r = zeros(1, Q);
    q(1) = sum(v(11:20)); % radius 10 .. 19
    r(1) = sum(t(11:20));
    for j = 2:Q
        q(j) = sum(v(11+5*j:15+5*j));
        r(j) = sum(t(11+5*j:15+5*j));
    end
    
    
    
    
    %% original
%    q(1) = sum(v(11:18)); % radius 10 .. 17
%    r(1) = sum(t(11:18));
%    for j = 2:Q % radius 18 .. 21 / ...
%        q(j) = sum(v(12+3*j:14+3*j));
%        r(j) = sum(t(12+3*j:14+3*j));
%    end
    q = q ./ r;
    % normalize
    % norm1 -> compare to maximum value
    q_norm1 = q / max(v);
    % norm2 -> sum is 1
    q_norm2 = q / sum(q);
    % output two sets
    ret1 = [ret1 q_norm1'];
    ret2 = [ret2 q_norm2'];
end

end